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[introduction 


Recent  seismological  investigations  of  nuclear  tests 


at  the  Nevada  Test  Site  (NTS)  have  manifested  the  need 


for  multidimensional  geophysical  models  of  the  shallow 


(less  than  5  km)  structure.  Careful  analysis  of  seismic 


signals  has  established  a  pattern  of  response  variations 


both  at  Yucca  Flat  and  Pahute  Mesa.  These  variations 


include  both  amplitude  and  travel  time  perturbations 


(Alewine,  personal  communication,  Minster  et  al.,  1981 


and  Taylor,  1983)  and  may  also  be  azimuth  dependent 


(Lynnes  and  Lay,  1984) . 


The  body  wave  magnitude  anomalies  observed  by 


Alewine  are  plotted  in  Figure  1  along  with  the  Bouguer 


gravity  anomaly.  The  gravity  anomaly  appears  to  reflect 


the  structure  of  the  14  my  old  Silent  Canyon  caldera 


collapse  feature,  which  is  now  buried  beneath  younger 


volcanics  (locally  to  over  2  km  depth)  .  The  body  wave 


magnitude  anomaly  correlates  rather  well  with  the  gravity 


feature,  which  strongly  suggests  a  common  explanation  for 


the  seismic  response  anomaly  in  terms  of  the  Silent 


Canyon  caldera  structure.  This  hypothesis  is  best 


approached  by  the  construction  of  a  structural  model 


using  geologic  and  geophysical  data  which  are  largely 


independent  of  the  nuclear  tests  themselves.  A  reliable 


model  of  this  type  can  be  used  to  test  the  seismic 


response  to  confirm  or  reject  the  hypothesis.  Ferguson 
(1987)  applies  this  approach  to  Yucca  Flat  with  some 
success.  It  is  important  to  explain  as  much  of  the 
variation  as  possible  with  known  (i.e.  independently 
determinable)  features  of  the  crust  and  then  relegate  the 
remaining,  reduced  variation  to  less  accessible 
explanations. 

Two-dimensional  structural  models  can  be  treated  by 
a  number  of  techniques,  such  as  finite  differences 
(Kelley  et  al.,  1976),  the  Aki-Larner  method  (Aki  and 
Larner,  1970)  and  generalized  ray  theory  (Helmberger  et 
al.,  1985)  to  name  a  few  examples.  These  methods  are 
being  applied  to  NTS  structures  (McLaughlin  et  al.,  1987 
and  Ferguson,  1987) ,  particularly  at  Yucca  Flat  where  the 
structure  is  fairly  simple  and  better  known  (Ferguson  et 
al.,  1988).  Techniques  applicable  to  three-dimensional 
models  may  be  practical  in  the  near  future. 

In  order  to  meet  the  new  requirements  for 
geophysical  modeling  at  Pahute  Mesa,  we  have  attempted  to 
gather  pertinent  geologic  and  geophysical  data  and 
produce  a  consistent  interpretation  of  these  data  in 
terms  of  geophysical  profiles.  The  use  of  wide  angle 
seismic  surveys  at  Pahute  Mesa  is  an  innovation  which  is 
being  tested  as  a  part  of  this  study.  The  approach  taken 
here  is  to  make  use  of  geologic  mapping  and  borehole 
explorations  to  a  full  extent.  Borehole  geophysics  is 


useful  in  the  calibration  of  geologic  models  for 
geophysical  purposes.  Density  and  velocity  measurements 
made  in  boreholes  have  been  assembled  and  analyzed. 
Gravity  observations  permit  a  uniform  areal  coverage 
where  boreholes  are  limited  and  are  often  a  good 
indication  of  seismically  important  structures.  Gravity 
data  are  availible  at  NTS  and  new  data  are  inexpensive  to 
obtain.  Seismic  reflection  and  refraction  have  not  been 
used  in  Pahute  Mesa  investigations  previously  due  to 
several  mitigating  factors.  A  novel  approach  to  the 
seismic  survey  has  been  attempted,  but  it  remains  to  be 
seen  whether  or  not  the  data  will  be  sufficiently 
sensitive  to  the  lateral  structure  variations  at  Pahute 


The  geology  of  the  Pahute  Mesa  has  been  the 
subject  of  numerous  publications  since  the  publication 
of  Nevada  Test  Site,  Memoir  110  of  the  Geological 
Society  of  America  (Eckel,  1968).  Important  papers  in 
that  volume  include  Noble  et  al.  (1968)  and  Orkild  et 
al.  (1968)  which  identify  and  describe  the  Silent 
Canyon  caldera.  Ten  years  of  geologic  mapping  on  7  1/2 
minute  quadrangles  was  complied  at  a  scale  of  1:48000 
by  Orkild  et  al.  (1969)  on  the  Geologic  Map  of  Pahute 
Mesa,  Nevada  Test  Site  and  Vicinity,  Nye  County, 
Nevada.  Radiometric  age  dates  were  published  by  Kistler 
(1968)  and  Marvin  et  al .  (1970).  Papers  on  the  Timber 
Mountain  and  related  calderas  to  the  south  were 
included  in  Memoir  110  and  were  also  published  in  the 
1970's  (Byers  et  al .  ,  1976  and  Christiansen  et  al., 
1977) .  Exploration  for  nuclear  waste  isolation  sites  to 
the  south  of  Timber  Mountain  and  continued  drilling  in 
support  of  the  testing  program  at  Pahute  Mesa  have 
encouraged  a  more  recent  synthesis  of  the  geologic  data 
by  Warren  et  al.  (1985).  Large  scale  geophysical 
structure  of  the  test  site  has  been  pursued  by  Minster 
et  al.  (1981),  Taylor  (1983)  and  Hoffman  and  Mooney, 
(1984).  This  sample  of  literature  is  not  exhaustive, 


but  is  sufficient  to  support  the  overview  of  the 
geologic  history  to  be  presented  here. 

The  southwestern  Nevada  volcanic  field  extends 
from  Goldfield,  Nevada  southeast  to  Mercury,  Nevada.  A 
number  of  major  volcanic  centers  occur  along  a  north- 
south  trend  along  the  western  boundary  of  NTS.  silent 
Canyon  caldera  is  north  of  the  most  persistent 
center(s)  at  Timber  Mountain  and  is  almost  totally 
buried  under  younger  volcanics  at  Pahute  Mesa  (Figure 
2)  . 

The  geologic  history  is  summarized  in  Table  1.  The 
volcanism  of  the  southwest  Nevada  volcanic  field  begins 
about  15  or  16  my  ago  with  eruptions  from  undefined 
centers  and  some  material  (the  Rhyolite  of  Kawich 
Valley)  from  Silent  Canyon  caldera.  These  older 
volcanics  overlie  Paleozoic  clastic  and/or  carbonate 
rocks  and  possible  granitic  rocks  similar  to  the  Climax 
and  Gold  Meadows  plutons  at  Pahute  Mesa.  The  volcanism 
may  be  related  to  the  occurence  of  north-south  trending 
normal  faults.  About  14  my  ago  the  Silent  Canyon 
caldera  collapsed  producing  the  Belted  Range  tuff. 
These  units  are  unusual  in  that  they  are  peralkaline  in 
composition  rather  than  the  more  common  calc-alkaline 
of  the  older  and  younger  volcanics.  We  hypothesize  that 
the  Belted  Range  tuff  and  stratigraphically  equivalent 
lavas  form  the  geophysical  basement  under  Pahute  Mesa. 


Figure  2.  The  major  volcanic  centers  of  the 
southwestern  Nevada  volcanic  field  are  shown  with  an 


Table  1 


Period 

IHHMII 

DescriDtior. 

L«C« 

Quaternary 

Alluvlua 

Ranging  froe  a  thin  veneer 
to  hundreds  of  asters  in 
the  grabene 

Basalt  lavas 

(Buckboard  Mesa) 

Locally,  saall  aeounts 

Early  Quaternary 
/Lac*  Tertiary 

Basin  and  range  graben 
development  on  MS  faults 
with  EH  extension  and  lata 
developing  strike-slip 
with  N30W  extension 

Lac*  Tertiary 

X-Ar  agea 

9. 3-6. 2  ay 

(Thirsty  Canyon) 

Black  Mt.  caldera 

li-lo. 4  ay 

(T labor  Ht) 

Timber  Mt.  caldera  collapse 

HZ  trending  faults  with  H78W 
extsnsion 

13-12  ay 

(Paintbrush) 

Claim  Canyon  caldaron 

MW  trending  faults 

13.2-12.4  ay 

calc-alkaline 
(Area  20) 

Sllant  canyon  caldera 
subaidenca,  no re  peralkaline 
on  the  eaec  side 

MS  faults  Inactive 

13.3-14  ay 

calc-aUcallna 
(Crater  Plat) 

crater  riat  caldera 

13.8  ay 

peralkaline 
(Belted  Range) 

Silent  canyon  caldera  collapse 

Morsel  faulting  on  MS  trend 

16-14  ay 

(Bedrock  Valley) 
(Fraction) 

(Tunnel  Beds) 
(Rhyolite  of 

Kavi.cn  valley) 

Local  volcanic  centers 

Silent  Canyon  caldera 

27  ay 

Inters ediate 
volcanic* 

(Monotony  Valley) 

Wide-spread  voicanisa  :n 
a  M£  trend  MTS  to  Utah 

Early  Tertiary 

High-angle  faulting 

MW  trend  (?) 

Lae*  Cr*tac*ous 

Granitic  intrusive* 
(Clisax  stock) 

Middle  Co  early 
Mesoroic 

Mo  sediaentation 

Orogeny 

Persian  to 
Pennsylvanian 

Carbonates 

Stable  aargm  (12  <n  total 
thickness  pc  to  now) 

Mississippian  to 
late  Devonian 

elastics  frea 
Roberts  Mt.  thrust 

Crogeny 

Early  Devonian  to 
late  Canbrian 

Carbonates 

Stable  targm 

Early  Caabrian  to 
laee  Precanrnan 

elastics 

St3bie  aargm 

*  X-Ar  ages  «roa  Kutier  (1968)  and  Marvin  at  al.  19*0) 


The  younger  Tuffs  and  Lavas  of  Area  2  0  were  erupted 
into  Silent  Canyon  caldera  about  13  my  ago  and  the 
caldera  was  filled  by  Paintbrush  tuff  from  Claim  Canyon 
cauldron,  on  the  south  side  of  Timber  Mountain,  between 
13  and  12  my  ago. 

The  collapse  of  Timber  Mountain  caldera  produced 
the  Timber  Mountain  tuff,  including  the  vitrophyre  of 
the  Rainer  Mesa  member,  which  caps  the  by  now  filled  in 
Silent  Canyon  structure.  Note  that  the  area  to  the 
southeast  of  Pahute  Mesa,  although  now  a  topographic 
low,  must  have  been  structurally  high  throughout  most 
of  this  time.  Belted  Range  rocks  are  now  exposed  there. 
The  exact  nature  of  the  structure  between  the  Silent 
Canyon  and  Timber  Mountain  calderas  is  also  not  well 
known.  A  small  amount  of  Thirsty  Canyon  tuff  from  the 
Black  Mountain  center  to  the  west  is  found  locally  at 
Pahute  Mesa.  Both  northwest  and  north-south  trending 
faults  have  been  active  at  Pahute  Mesa  during  the 
Tertiary  and  northeast  trending  faults  have  been 
important  regionally.  During  the  Paintbrush  deposition, 
the  northwest  trending  Silent  Canyon  shear  zone  was 
active  in  northern  Pahute  Mesa,  truncating  the  north- 
south  faults  and  bounding  the  low  lying  depositional 
area.  The  Split  Ridge,  Scrugham  Peak,  Almendro,  West 
Greeley,  East  and  West  Boxcar  and  Purse  north-south 


trending  faults  have  been  significant  throughout  the 
caldera  history  (Figure  3) . 

The  caldera  collapse  entailed  a  more  or  less 
elliptical  zone  of  faulting  or  slumping  as  indicated  on 
the  map  of  Orkild  et  al.  (1969)  and  in  Figure  3.  The 
margin  is  not  well  exposed  and  the  boundary  shown  is 
interpretive  in  nature.  Geophysical  investigations  may 
help  to  clarify  the  margin  structure  and  the 
configuration  of  the  post-Belted  Range  Silent  Canyon 


caldera  structure. 


Borehole  geophysical  logs  are  routinely  made  at  NTS 
in  exploratory  and  emplacement  holes  for  containment 
studies.  In  this  investigation  borehole  gravity,  density 
logs  and  borehole  velocity  surveys  have  been  of 
particular  interest.  Well  logs  from  23  holes  were 
obtained  from  the  GEODES  data  base  at  Los  Alamos.  The  log 
availablity  is  summarized  in  Table  2.  Only  two  wells  in 
the  western  half  of  the  Mesa  were  studied  and  the  depth 
extent  is  generally  limited  to  less  than  700  m.  Where 
deep  (1  to  2  km)  holes  exist,  only  density  information  is 
available.  The  important  holes  used  in  this  investigation 
are  plotted  on  Figure  4. 

Several  goals  were  identified:  1)  Can  geophysically 
distinct  geologic  units  be  identified  ?  2)  Can  the 
geologic/geophysical  units  be  characterized  by  simple, 
well  defined  velocity  and  density  functions  ?  3)  Can  a 
relationship  be  found  between  velocity  and  density  so 
that  velocity  can  be  predicted  from  the  more  abundant 
density  data  ? 

It  should  be  pointed  out  that  this  study  has  been 
hampered  by  a  lack  of  good  stratigraphic  logs  for  these 
holes.  Some  of  the  older  logs  are  in  error  and  positive 
identification  of  some  of  the  units  is  difficult  at  best. 
Recent  studies  by  Richard  Warren  at  Los  Alamos  using 


Well 

Density 

Gravity 

Velocity 

Water 

Level 

Depth 

Extent 

U19aa 

* 

* 

709 

600 

U19af 

* 

* 

* 

695 

Ul9ae 

* 

* 

* 

707 

700 

U19ak 

* 

* 

655 

600 

Ul9ab 

★ 

* 

* 

700 

Ul9x 

* 

* 

* 

715 

700 

Uel9z 

* 

* 

668 

800 

U20n 

* 

625 

1200 

PM-1 

* 

644 

2200 

Ue20f 

* 

1300 

U19ad 

* 

700 

Uel9n 

* 

400 

U19x 

* 

700 

Uel9pl 

* 

i 

700 

U19g 

* 

* 

627 

700-100 

U19ar 

* 

700 

Uel9p 

* 

690 

600 

Uel9gs 

* 

1300-23 

Uel9fsl 

* 

800-140 

Uel9i 

* 

900 

Uel9x 

* 

800 

Uel9g 

* 

i 

627 

700-900 

indicates  that  the  log  is  available. 


All  depths  are  in  meters 


geochemical  and  petrologic  techniques  have  helped  to 
eleviate  this  problem,  but  the  work  is  as  yet  incomplete. 

Initially  the  borehole  gravity  in  the  deepest  holes, 
like  PM-1  and  Ue20f,  were  analyzed  to  obtain  a  density 
function  for  the  caldera  filling  tuff  units.  The  function 

a(z)  =  1750  +  420  z,  (km,  kg/m3) 

was  fit  to  PM-1.  Exceptions  to  this  relation  are 
higher  density  lava  units  and  the  Belted  Range  tuff  (Tbt) 
which  is  associated  with  what  we  term  the  pre-caldera 
basement.  The  Tbt  appears  to  have  a  density  near  2  500 
kg/m3.  This  model  was  found  to  predict  the  tuff  densities 
below  0.5  or  0.6  km  in  other  borehole  gravity  logs.  The 
model  and  the  density  data  for  PM-1  and  Ue2  0f  are  shown 
in  Figure  5.  Borehole  gravity  is  a  very  good  way  to 
characterize  bulk  density  due  to  its  sensitivity  to  a 
volume  of  rock  around  the  borehole,  for  the  same  reason 
it  can  also  be  subject  to  bias  if  significant  structures 
are  nearby  (Hearst  and  Carlson,  1982). 

The  holes  containing  both  density  and  velocity  data 
were  analyzed  in  greater  detail,  particularly  Ul9af, 
U19ae,  U19x  and  U19ab  which  are  near  the  seismic  line. 
The  more  distant  holes  Uel9z  and  ul9ak  were  also  studied. 
Geophysical  data  for  these  holes  were  grouped  according 
to  rock  type  by  intervals.  These  groups  were  Rainer  Mesa 
tuff  (Tmr) ,  all  other  tuffs  and  lavas.  Zeroth  and  first 
degree  polynomial  models  were  fit  to  the  data  for  each 


B 

&888mm 

55 

Figure  5.  Densities  computed  from  borehole  gravity  in 
wells  PM-1  and  Ue20f.  The  density  function  derived  from 
PM-1  below  600  m  is  shown  along  with  the  density 
profiles. 


krtarval  Danatty  fro/m3) 


hole  by  group  and  for  all  holes  combined  by  group  by 
least  squares.  It  should  be  noted  that  the  use  of  least 
squares  for  this  purpose  is  probably  not  appropriate,  but 
the  results  seem  to  be  consistent.  The  models  based  on 
the  data  combined  from  all  six  holes  are  summarized 
below: 


Velocity  (km  and  km/sec) : 

Tuffs  (including  Tmr)  - 
v ( z )  =  1.70  +  2.15  z 
Lavas  - 
v ( z )  =  2.9 

Density  (km  and  kg/m3)  : 

Tuffs  (excluding  Tmr)  - 
a (z)  =  1510  +  740  z 
Lavas  and  Tmr  - 
o( z)  =  2200 

The  root  mean  square  error  for  the  velocity  models 
is  0.7  to  1.0  km/sec  and  100  to  120  kg/m3  for  the  density 
models.  The  correlation  coefficient  is  an  acceptable  0.8 
for  the  tuff  density  model  and  a  rather  poor  0.5  for  the 
tuff  velocity  model.  The  models  are  plotted  along  with 
logs  for  U19af  and  Uel9z  in  Figure  6.  The  lavas  are  not  a 
well  defined  geophysical  unit.  The  lava  intervals  are 
somewhat  indistinct  and  the  physical  properties  are 
highly  variable.  They  do  show  a  tendency  to  a  higher 
density,  but  the  velocity  is  not  well  distingished  from 


Figure  6.  Density  and  P-wave  velocity  profiles  from 
U19af  and  Uel9z.  The  shallow  velocity  and  density 
models  are  also  plotted. 


Velocity: 

Tuffs- 

v(z)  =  1.70  +  2.15  z 
Lavas- 

v(z)  =  3.50  km/sac 
Density: 

Tuffs  (excluding  Tmr)- 

p(z)  =  1.51  +  0.74  z 


km/sec 

z  in  km 


gm/cm3 


Lavas  and  Tmr- 


p(z)  =  2.20  gm/cm* 
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the  tuffs.  The  lava  velocity  distribution  is  better 
characterized  by  a  large  variance  and  a  broad,  almost 
uniform  probability  distribution.  If  the  lavas  were  to  be 
grouped  with  the  tuffs  the  result  would  be  effected  only 
slightly.  The  Rainer  Mesa  tuff  is  similar  in  velocity 
structure  to  the  other  tuffs  and  is  definitely  of  a 
higher  density. 

The  velocity  and  density  models  developed  for  these 
six  holes  are  all  based  on  data  acquired  above  the  water 
table,  which  occurs  at  about  0.65  to  0.7  km.  The  density 
model  developed  for  PM-1  is  valid  only  below  depths  of 
0.6  km.  It  is  significant  that  the  intercept  density  is 
greater  and  the  gradient  is  lower  for  the  deep  model. 
Some  implications  of  this  will  be  discussed  in  a  later 
paragraph.  A  similar  effect  should  also  occur  in  the 
velocity  vs.  depth  function.  A  velocity  model  recently 
published  by  Leonard  and  Johnson  (1987)  supports  this 
conclusion.  Their  model  is  displayed  in  Figure  7  along 
with  the  tuff  velocity  model  developed  here  and  the  raw 
interval  velocities  for  tuff  and  lava  intervals.  The 
linear  polynomial  model  is  biased  slightly  low,  but  the 
gradient  is  similar  down  to  0.6  km  depth.  The  Leonard  and 
Johnson  model  was  found  by  a  linearized  inversion  of  the 
travel  times  from  nuclear  tests  out  to  offsets  of  12  km. 
The  model  represents  a  gross  vertical  and  horizontal 
average  of  the  true  velocity  distribution.  The  tremendous 


Figure  7.  The  shallow  velocity  model  from  this  report, 
the  Leonard  and  Johnson  velocity  profile  and  the  raw 
interval  velocities  from  Area  19  holes  are  plotted 
together  on  this  figure  for  comparison. 


scatter  of  the  the  measured  interval  velocities  about 
these  average  models  should  be  noted.  At  a  depth  of  2  km 
the  velocity  approaches  4.6  km/sec  and  the  travel  time 
data  lose  resolution.  A  clear  break  in  the  velocity 
gradient  is  observed  in  the  vicinity  of  the  water  table. 
The  deeper  velocities  are  predicted  by 
v(z)  =  3.1  +  0.6  z  km/sec. 

Geophysically,  zones  of  wet  and  dry  tuff  are  distingished 
on  the  basis  of  property  gradient  with  depth. 

The  density  and  also  the  velocity  of  clastic 
sedimentary  rocks  depend  greatly  on  the  porosity  and 
these  volcanic  rocks  might  also  be  expected  to  display 
this  effect.  The  tuffs  are  initially  very  porous  (30  %  or 
more)  and  the  porosity  decreases  with  depth  due  to 
compaction.  The  compaction  and  porosity  are  also  strongly 
effected  by  secondary  mineralization,  glassiness  of  the 
matrix  and  welding.  The  lavas  could  have  a  great  range 
of  properties  from  vugginess  to  very  low  porosity.  We  can 
make  some  speculations  based  on  a  theoretical  compaction 
model  for  sedimentary  rocks  (Korvin,  1984)  .  The  porosity 
should  obey  an  exponential  law  of  the  form 
0  =  0O  exp(-kz)  . 

If  k  is  small  this  might  be  reasonably  replaced  by  a 
linear  relationship,  such  as 
0  =  0O  -  kz  . 

Given  a  simple  mixing  law  for  the  density 
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we  can  obtain  the  linear  density  function 

c(z)  =  am  -<Pq  (cjjjj  -  of)  +  k(om  -  af)z  . 

If  we  assume  that  the  porosity  decrease  rate,  k,  and  the 
matrix  density,  om,  are  constant  and  that  the  only 
difference  in  the  shallow  and  deep  density  functions  is 
the  pore  fluid  density.  Of,  we  can  conclude  that  k  =  0.28 
or  0.29,  for  air  in  the  shallow  pores  and  water  in  the 
deeper  pores.  The  difference  in  the  intercept  values  for 
the  density  models  is  less  well  explained  by  this 
hypothesis,  as  a  discrepancy  of  160  kg/m3  occurs.  The 
surface  porosity,  0O,  is  found  to  be  4  0  %,  which  is 

reasonable  for  these  tuffs. 

The  velocity  models  are  much  more  difficult  to 
reconcile  with  the  porosity.  Published  empirical 
relationships  (Gardner,  et  al.,  1974)  between  porosity 
and  velocity, 

v  =  <f>  /  vf  +  (1-0)  /  vm  , 
and  density  and  velocity, 

<j  =  a  v^ *25  ^ 

are  inconsistent  and  neither  seems  to  fit  the  density  or 
porosity  functions  estimated  from  the  density  data.  It 
was  expected  that  the  shallow  and  deep  velocity  functions 
could  be  explained  by  a  simple  change  in  pore  fluid 
similar  to  the  density  relationships,  but  this  has  not 
been  realized.  The  empirical  relationship  in  Ibrahim  and 


Keller  (1981)  for  volcanic  rocks  does  not  fit  these  data 
and  it  is  difficult  to  tell  if  the  rocks  studied  in  that 
paper  are  at  all  similar  to  those  at  Pahute  Mesa.  Due  to 
this  problem  a  predictive  relationship  between  density 
and  velocity  has  not  been  established  and  a  different 
approach  will  be  taken  in  the  future.  The  availability  of 
the  Leonard  and  Johnson  velocity  model  makes  the 
prediction  of  velocity  from  density  less  crucial  to  the 
geophysical  modeling  effort. 


Gravity  Modeling  Technique 
Fourier  Gravity  Inversion  Technique: 

The  gravity  modeling  technique  to  be  applied  at 
Pahute  Mesa  must  be  capable  of  using  the  linear  density 
function  indicated  by  the  log  analysis.  Although 
initially  two-dimensional  models  will  be  considered, 
Three-dimensional  calculations  should  be  possible. 
Ferguson  et  al .  (1988)  have  demonstrated  that  the  so 
called  Parker-Oldenburg  technique  (Parker,  1972 
Oldenburg,  1974)  is  effective  in  the  inversion  of  large, 
three-dimensional  problems.  In  this  report  an  extension 
of  the  original  formulation  to  include  a  linear  density 
function  will  be  made.  Although  no  three-dimensional 
calculations  will  be  made  here,  the  extension  from  two  to 
three  dimensions  is  not  difficult. 

Parker  (1972)  formulated  the  forward  potential 
field  calculation  by  assuming  the  source  of  the 
gravitational  anomaly  to  be  a  layer  at  some  depth,  z0, 
bounded  by  variable  interface  topography  h(r)  and  g(r), 
and  density  contrast  o(£)  (Figure  8).  Depth,  z,  is  taken 
to  be  positive  downward.  Starting  with  the  Fourier 
transform  of  the  gravitational  potential  given  by  : 


F[U(r0)  =  G  o 


.h(r) 


exp(i  k-r)  /  | r 0  —  r I  dz  dS  ds( 


D  g  (£) 


rQ  =  position  of  the  observation 


Figure  8.  Schematic  description  of  parameterization  of 
the  Fourier  gravity  inversion  technique  as  derived  by 
Parker  (1972)  and  Oldenburg  (1974).  Note  that  in  this 
description,  z  is  positive  downward  and  o  is  a  function 
of  lateral  position. 


r  =  position  within  the  volume  to  be  integrated 
£  =  projection  of  r  onto  the  surface  z  =  0 
and  integrating  over  domains  X  and  z,  the  resulting 
integral  becomes: 


F[U(r0)]  =  2f  G  a 


exp(i  k-2L  -  Ulz0) 


♦  {  exp(|)s|  (h(£)  -  g(£)))  }  /  |k|2  dS 


After  a  Taylor  series  expansion  of  the  exponential 
function  of  the  interface  topography,  exp[h(r) -g(r) ] , 

00 

F[U(r0)]  «  2ir  G  a  exp(-|&|z0)  •  £  lKln"2  /  n! 

n=l 

*  F[hn(£)  -  gn(£)] 


the  integral  then  becomes  a  Fourier  transform  of  sums  of 
successive  powers  of  the  interface  topography. 
Differentiating  this  expression  yields  the  equation  for 
the  gravity  effect: 

00 

F[6g(£)  ]  =  2ir  G  o  exp(-|fc|z0)  *  £  l&ln_1  /  n! 

n=l 

•  F[hn(£)  -  gn(£) } 

Oldenburg  (1974)  derived  the  inverse  solution  by 
solving  for  the  first  order  source  topography  term, 
[h (£) -g (£) ] .  This  results  in  the  inverse  solution  by  the 
method  of  successive  approximations: 


F[h(£)  -  g{£) ]  =  -F[Sg(c)J  exp(|k|z0)  /  (2n  G  o) 
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ao 

•  2  |jj|n-1  /  n!  F[hn(£)  -  gn(r)  ] 
n=2 

The  derivation  by  Parker  (1972)  and  Oldenburg  (1974) 
could  be  formulated  with  a  laterally  varying  density 
function,  a(i) : 


F[a(£)  (h(£)  -  g(r))]  =  -F[Sg(r)]  exp(|k|z0)  /  (2tt  G  a) 

00 

•  2  I&!n_1  /  n!  F[0(r)  (hn(r)  -  gn(r))]  . 

n=2 

This  result  is  a  non-linear  inverse  solution  for 
parameterization  of  a  depth  model  using  potential  field 
data.  As  discussed  by  Oldenburg  (1974) ,  parameterization 


of  this  model  is  non-unique.  h  priori  specification  of 


the  density  distribution  and  one  of  the  bounding  surfaces 


reduces  some  ambiguity  of  this  non-uniqueness. 


Addition  of  a  Vertically  Varying  Density  Function: 


Since  density  is  a  function  of  depth  within  the 


Earth,  it  is  useful  to  derive  the  Fourier  transform  of 


the  gravity  effect  with  a  vertically  varying  density 


function.  Assume  that  c(z)  is  a  linear  function  such 


that 


a(z)  »  an  +  cz  , 


where  o0  is  the  density  contrast  on  the  observation 


plane  and  c  is  the  density  gradient.  Replacing  o  in  the 


forward  calculation  with  o(z)  results  in  the  following 


expression  for  the  two-dimensional  case, 


\ 4 


— b 


F[Sg(£)]  =  21T  G  exp(-|fc|z0)  {  oQ  £  (k|n_2  /  n! 

n=l 

00 

•  F[hn(r)  -  gn(£)  ]  +  C  E  | is | n_1  /  n! 

n=2 

•  F[hn(£)  -  gn(E) ] 

where  z  is  taken  to  be  positive  downward.  The  inverse 
expression  then  becomes: 

F[h(r)  -  g  (£)  ]  =  F[6g(r)]  exp(|k|z0)  /  (2tt  G  o) 

00 

-  2  {  lfcln_1  /  n!  +  (c/o0)  \K\n~2  /  n!  (n-1)  } 
n=2 

•  F[hn(r)  -  gn(r) ] 

A  detailed  derivation  of  this  result  is  given  in  the 
Appendix.  Note  that  instead  of  summation  over  one 
infinite  series,  this  result  contains  two  infinite 
series.  The  first  series  is  essentially  the  same  as  that 
in  the  original  results  of  Parker  (1972)  and  Oldenburg 
(1974).  The  second  series,  which  adds  the  effects  of  the 
vertical  density  gradient,  has  no  apparent  effect  on  the 
first  order  interface  term,  which  is  the  largest  term, 
since  series  summation  begins  at  n  =  2.  The  effect  of  a 
density  gradient  on  the  first  infinite  series  becomes 
apparent  by  considering  the  effect  of  shifting  the 
observation  plane,  zQ  ,  as  suggested  by  Parker  (1972)  to 
speed  convergence  of  the  solution.  For  the  constant 
density  case  this  shifting  does  not  affect  the  magnitude 


of  the  resulting  gravitational  anomaly.  In  the  case  of 
the  vertical  density  function,  if  the  zQ  plane  is 
shifted,  the  initial  density,  aQ,  must  also  be  shifted. 
Magnitude  of  the  gravity  effect  in  the  forward  solution 
and  model  depth  in  the  inverse  solution  are  then  scaled 
as  a  result  of  the  density  gradient. 

Since  each  series  is  summed  separately,  the 
convergence  of  each  series  should  be  tested  separately. 
Tests  for  convergence  of  the  Taylor  series,  as  discussed 
by  Parker  (1972)  are  applied.  Convergence  of  the  inverse 
solution  is  also  controlled  by  the  stability  of  the 
explicit  downward  continuation  function,  exp(|&|z0). 

Shallow  source  regions  and  observational  error 
cause  the  continuation  function  to  increase  rapidly  with 
increasing  depth  and  force  the  inverse  solution  to 
diverge.  Low  pass  filtering  of  the  data  suppresses  this 
effect.  However,  eliminating  short  wavelength  effects 
with  an  arbitrary  low  pass  filter  results  in 
overfiltering  or  underfiltering  of  the  data  introducing 
greater  misfit  between  observed  and  theoretical  gravity 


data 

(Ferguson 

et  al .  , 

1988) . 

In 

addition,  if 

the 

bandwidth  of  the 

filter 

is  too 

large , 

the 

inclusion 

of 

high 

wavenumber 

information 

will 

still 

result 

in 

divergence  of  the  inverse  solution. 

A  better  method  of  stabilization  of  the  downward 
continuation  process  can  be  achieved  by  application  of  a 
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regularization  function  as  described  by  Tikhonov  et  al. 
(1968)  : 

f(&,a,z)  =  1  /  (1  +  a2  ]c2  exp(lKlz)) 

The  effect  of  this  filter  is  to  taper  the  downward 
continuation  function  before  it  begins  to  grow  large. 
For  low  wavenumber  values,  the  filter  has  a  minimal 
effect  (Figure  9) .  As  the  wavenumber  increases,  the 
exponential  function  begins  to  decay  rather  than  grow. 
The  amount  of  tapering  is  controlled  by  the  parameter  a. 
As  shown  in  Figure  9,  application  of  the  filter  with  a 
large  a  value  rapidly  tapers  the  exponential  function 
whereas  a  more  gradual  taper  is  achieved  with  smaller  a  . 
Optimizing  the  regularization  filter  consists  of  finding 
an  optimal  value  for  a  which  controls  the  growth  of  the 
downward  continuation  without  tapering  too  rapidly. 
Ferguson  et  al.  (1988)  used  the  Nelder-Mead  polytope 
technique  (1965)  to  search  the  solution  space  for  the  a 
value  which  best  minimized  the  root  mean  square  (rms) 
error  between  observed  and  calculated  gravity  anomaly 


values  after 

each 

iteration 

of  the  inversion. 

A 

limitation  of 

this 

method  is 

that  a 

reasonable 

first 

approximation 

for 

a  is  required 

to  insure 

rapid 

convergence  of 

the 

Nelder-Mead 

simplex 

technique . 

Also , 

computational  efficiency  of  the  Fourier  gravity  inversion 
technique  is  decreased  because  of  additional  evaluations 


Figure  9.  Amplitude  of  the  downward  continuation 
function  increases  exponentially  with  increasing 
wavenumber.  Regularization  tapers  the  exponential 
growth  with  the  amount  of  tapering  proportional  to  the 
value  of  the  filter  parameter  a.  This  figure 
illustrates  the  effect  of  increasing  a  values  on 
downward  continuation. 
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of  the  forward  solution  used  in  the  search  for  an  optimal 
value  for  a. 

Characteristics  of  Potential  Field  Spectra: 

To  optimize  of  the  regularization  filter,  spectral 
characteristics  of  potential  field  data  are  considered. 
Spector  and  Grant  (1970)  and  Hahn  et  al.  (1975) 
demonstrated  that  potential  field  spectra  have  a 
characteristic  shape  which  is  dominated  by  the  effect  of 
source  depth.  The  effect  of  source  depth  produces  a 
linear  relationship  between  the  logarithmic  power 
spectrum  and  wavenumber,  the  slope  of  the  line  being 
proportional  to  the  source  depth  (Figure  10)  .  A  power 
spectrum  of  gravity  anomaly  data  produced  by  sources  at 
multiple  depths  can  be  characterized  by  a  series  of  line 
segments  with  different  slopes. 

Shallower  sources  have  correspondingly  lower 
spectral  slopes.  Observational  errors,  if  statistically 
uncorrelated  to  the  data,  are  represented  by  "white 
noise,"  the  mean  spectral  value  of  which  is  a  constant  or 
flat  line.  As  is  shown  in  Figure  10,  the  effects  of 
deeper  sources  dominate  the  low  wavenumber  portion  of  the 
spectrum,  whereas  shallower  sources  and  white  noise 
dominate  the  higher  wavenumbers.  In  this  study,  the 
gravity  effect  of  deeper  sources  is  considered  to  be  the 
anomaly  of  interest;  therefore,  the  portion  of  the 
spectrum  represented  by  deeper  sources  is  assumed  to  be 


mmmz- 


Figure  10.  Graphical  representation  of  potiental  field 
spectra.  Part  (a)  illustrates  a  spectrum  with  sources 
at  three  different  depths,  Zj_,  z2,  and  z3.  Slopes  of 
line  segements  representing  each  depth  are  proportional 
to  the  depth  with  steepness  of  slope  indicating  greater 
depth.  Noise  is  indicated  by  the  flat  part  of  the 
spectrum.  Part  (b)  illustrates  the  assumption  used  in 
the  regularization  process;  a  single  source  at  depth  z 
produces  the  observed  gravity  effects  and  is  considered 
the  signal  (dominating  the  low  wavenumbers)  and  is 
distinct  from  noise  which  has  a  flat  spectrum 
(dominates  higher  wavenumbers) .  Crossover  wavenumber 
separates  signal  from  noise  and  occurs  when  the  ratio 


the  signal.  Other  portions  of  the  spectrum  representing 
shallow  sources  or  fluctuations  due  to  observational 
error  are  considered  undesirable  components  of  the 
overall  gravity  effect  and  are  assumed  to  be  noise. 

Since  amplitude  or  power  of  the  signal  decays 
linearly  for  potential  field  data  plotted  on  a  semi- 
logarithmic  scale,  the  effect  of  noise  begins  to  dominate 
as  the  slope  of  the  signal  falls  below  the  noise  level. 
The  point  at  which  the  signal/noise  ratio  equals  unity, 
termed  the  crossover  wavenumber  (kc) ,  is  the  wavenumber 
after  which  the  effect  of  noise  dominates  the  spectrum 
(Figure  10) . 

In  order  to  regularize  the  downward  continuation 
function  without  overfiltering  or  underfiltering  the 
data,  is  necessary  to  choose  the  regularization  parameter 
a  which  tapers  the  exponential  function  at  the  crossover 
wavenumber,  therefore  retaining  the  maximum  amount  of 
signal  with  a  minimum  amount  of  noise. 

Least  Squares  Spectral  Model: 

The  power  spectral  estimates  for  this  study  were 
made  by  calculating  the  periodogram  of  the  data  defined 
by: 

N-l 

IN(k)  =  ( 2/N)  |  z  6g(n)  exp(-ikn)  |2 

n=0 

where  k  *  wavenumber  (radians/km) 

N  *  number  of  observations 


Sg(n)  =  gravity  observations  (raGals) 

n  =  distance  index. 

Sx  *  the  spatial  sample  interval 

and  x(n)  =  n  •  6x,  where  n  =  0,1,..., N-l. 

The  periodogram  is  normalized  by  N/2  and  sampled 
at  a  discrete  set  of  wavenumbers,  kj ,  determined  by  the 
sampling  rate  of  the  data  as  follows: 

Sk  =  2ir  /(N  Sx)  , 

where  Sk  =  wavenumber  sampling  interval 

and  kj  =  j  •  Sk  where  j  =  0,1,..., N/2. 

In  all  periodograms  calculated,  the  data  are  scaled  so 
that  the  wavenumber  at  the  Nyquist  *  ir . 

The  periodogram  is  a  natural  but  unsatisfactory 
estimate  of  the  power  spectrum  due  to  its  fluctuating 
form  (Priestly,  1981) .  However,  in  this  study  the  power 
spectrum  is  reduced  to  a  few  rough  parameters,  such  as 
crossover  wavenumber  and  low  wavenumber  intercept,  and 
the  periodogram  is  found  to  be  an  adequate  spectral 
estimate . 

In  order  to  find  the  crossover  wavenumber,  a 
technique  for  distinguishing  between  signal  and  noise  is 
needed.  Since  the  spectral  shape  is  composed  of  a 
sloping  line  representing  the  signal  or  source  and  a  flat 
line  representing  noise,  a  least  squares,  line  spline 
approximation  to  the  spectrum  seemed  appropriate. 


The  signal  portion  of  the  spectrum  is  fitted 


separately  from  the  noise.  The  noise  is  constrained  to 


have  zero  slope,  and  both  solutions  are  constrained  to 


merge  at  the  crossover  wavenumber.  The  correct  choice 


for  kc  produces  the  best  least  squares  fit  to  the 


spectrum.  Optimization  for  kc  using  this  least  squares 


fitting  is  accomplished  with  a  univariate  minimization 


procedure,  FMIN  (Forsythe  et  al.,  1977). 


In  practice,  it  is  necessary  to  weight  the  least 


squares  fit  to  the  signal.  A  weighting  scheme  with  the 


spectral  signal  weight  being  inversely  proportional  to 


wavenumber  and  decreasing  linearly  works  well. 


weighting  factors  are  determined  a 


X  with  the  weight 


at  k  =  o  equal  to  one  and  the  weight  at  the  Nyquist  is 


equal  to  zero.  Since  the  amount  of  noise  in  different 


data  sets  varies,  the  weighting  factors  may  need  to  be 


adjusted  accordingly. 


Wavenumber  Dependence  of  the  a  Parameter: 


Once  a  method  for  determining  kc  is  found,  a 


method  for  choosing  the  regularization  filter  with 


rolloff  at  the  crossover  wavenumber  depends  on  finding  a 


relationship  between  the  regularizing  parameter,  a,  and 


kc.  Some  preliminary  investigation  of  the  behavior  of 


the  parameter  a  shows  that  for  regularization  at  a  given 


depth,  increasing  a  by  factors  of  10  decreases  the  rate 


of  the  filter  rolloff  (over  a  range  of  10  dB)  by  a  factor 


of  two  for  a  given  depth  (Figure  11) .  This  rough 
correlation  of  a  and  wavenumber  is  not  observed  for 
either  relatively  small  or  large  values  of  a.  In  the 
case  of  small  a,  rolloff  is  not  observed;  rather,  the 
spectrum  is  flat.  The  range  of  a  which  produces  flat 
spectra  decreases  with  increasing  depth  as  shown  by  plots 
of  filter  width  as  a  function  of  a  for  different  depths 
(Figure  12) .  For  large  a  values,  rolloff  occurs  at  very 
low  wavenumbers.  Therefore,  it  appears  that  for  a  given 
depth,  a  is  a  function  of  wavenumber,  and  a  relationship 
to  kc  can  be  established. 

Empirical  Relationship  of  a  to  kc  as  a  Function  of 
Continuation  Depth: 

A  procedure  is  developed  to  produce  statistical 
data  which  will  establish  the  empirical  relationship 
between  a  and  kc.  A  series  of  synthetic  gravity  models 
are  generated  at  several  depths  and  the  gravity  effects 
of  the  models  calculated.  The  model  used  initially  is  a 
line  source  since  the  gravity  effect  is  known  and  easy  to 
calculate.  In  addition,  the  spectrum  is  not  distorted 
by  the  effects  of  source  width  and  depth  extent  (Spector 
and  Grant,  1970) . 

The  gravity  effect  of  the  synthetic  model  is 
calculated  at  depths  of  z  =  0  (the  observation  plane)  and 
z  =  dc,  where  dc  is  the  continuation  depth,  taken  to  be 
0.9  times  source  depth.  Gaussian  random  noise  is  then 
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Figure  11.  Plot  of  log10  of  amplitude  of  the 
regularization  filter  as  a  function  of  wavenumber. 
Each  plot  from  0  to  9  represents  a  different  a  value, 
a  values  are  incremented  by  a  factor  of  10. 
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Figure  12.  Plot  of  regularization  filter  width  as  a 
function  of  the  log1Q  of  a  for  four  different  depths 
represented  by  plots  (a)  through  (d)  .  Filter  width  is 
defined  as  the  wavenumber  at  which  amplitude  of  the 
filter  decays  by  1/e. 
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added  to  the  surface  data,  i.e.,  superimposed  on  the  data 
are  an  additional  percentage  of  normally  distributed 
random  numbers.  The  maximum  gravity  anomaly  value,  here 
defined  as  the  signal  level,  is  used  to  scale  the  noise. 
Note  that  noise/signal  ratio  in  this  case  is  defined  as 
the  ratio  of  standard  deviation  of  the  Gaussian  random 
numbers  to  the  maximum  anomaly  value  and  is  applicable  in 
in  the  spatial  domain.  A  distinction  is  therefore  made 
between  spatial  noise/signal  ratio  and  noise/signal  ratio 
defining  crossover  wavenumber  in  the  wavenumber  domain. 
The  procedure  FMIN  is  again  used  to  optimize  the  choice 
for  a  by  minimizing  the  least  squares  error  between  the 
surface  data  after  regularization  and  downward 
continuation  and  the  actual  gravity  values  calculated  at 
the  continuation  depth.  The  crossover  wavenumber  is 
chosen  by  the  spectral  modeling  procedure  described 
previously. 

Depths  of  synthetic  models  range  from  1.25  km  to 
18  km?  depth  of  the  source  is  doubled  after  each  set  of 
calculations.  At  each  depth,  four  different  noise  levels 
are  applied  to  the  data;  noise/signal  ratios  are  0.02, 
0.05,  0.10,  and  0.20.  Figure  13  is  an  algorithm 
outlining  the  experimental  procedure.  Figure  14 
illustrates  the  results  of  execution  of  this  procedure 
for  a  line  source  at  6  km  (Figure  15)  and  a  noise/signal 
ratio  of  5  percent.  Note  that  the  anomaly  calculated  at 


‘For  each  j,  the  depth  of  the  source  is  doubled;  the  range  of  depths  is  from 
0.125  to  18  km* 

Generate  synthetic  gravity  data  at  the  surface  and  at  the  continuation  depth 
(0.9  times  the  model  depth)  for  a  line  source  at  depth  =  z-. 

For  i  =  1  to  4  do 

noise/signai  (n/s)  ratio: 
if  i  =  1  then  n/s  =  2  percent 

if  i  =  2  then  n/s  =  5  percent 

if  i  =  3  then  n/s  =  10  percent 

if  i  =  4  then  n/s  =  20  percent 

•noise  to  signal  ratio  =  ratio  of  the  standard  deviation  of  the  Gaussian 
random  noise  to  the  maximum  gravity  anomaiy  vaiue* 

Add  noise  to  the  surface  data 

Calculate  power  spectral  estimate  of  surface  data 

Optimize  for  correct  a  value  by  minimizing  least  souares  error 
between  regularized  surface  data  and  data  at  the  continuation  oectn 
using  a  nonlinear  optimization  procedure. 

Find  kc  by  least  squares  spectral  modeling. 

Filter  the  data  using  the  regularization  and  downward  continuation 
functions. 

Calculate  leas.  s  error  between  the  filtered  surface  lata  anc  :ata 

at  the  continue  .pth. 

Plot  results. 

End  for 

End  for 

Perform  linear  regression  analysis  on  data  from  ail  deptns 
Plot  a  as  a  logarithmic  function  of  l<c 

Plot  the  siope  of  the  lcg10  of  a  as  a  function  or  <c  .nes  as  a  ;  inct.cn  or 

depth  and  perform  .inear  regression 

Plot  the  Oq  intercept  of  the  log10  of  a  as  a  function  of  <c  .:nes  as  a 
function  of  depth  and  perform  linear  regression 

Derive  empirical  relationship  between  a  and  <c  as  a  net, on  of  cectn 


Figure  13.  Experimental  procedure  to  find  relationship 
between  the  regularization  parameter  a  and  crossover 
wavenumber  kc  as  a  function  of  depth. 
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Figure  14.  An  example  of  results  of  the  experimental 


procedure  described  in  Figure  13  for  a  line  source  at  a 


single  depth:  (a)  surface  data  with  5  percent 


noise/signal  ratio,  (b)  regularization  filter  with 


optimal  a  value  of  0.182,  (c)  comparison  of  regularized 


downward  continued  data  with  data  calculated  at  the 


continuation  depth  (0.9  times  the  source  depth)  with 


rms  error  of  0.0683  mGals,  9d)  least  squares  spectral 


model  with  crossover  wavenumber  at  0.688  radians/km. 
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the  continuation  depth  closely  approximates  a  delta 
function,  which  is  expected  for  a  line  source.  Even 
though  128  samples  were  used  to  define  the  anomaly,  less 
than  5  percent  of  the  total  number  of  samples  defined  the 
peak  amplitude.  Therefore,  even  though  the  filtered 
signal  fits  the  "true"  signal  at  the  continuation  depth 
in  a  least  squares  sense,  there  is  a  large  maximum  error 
between  peak  amplitudes  as  the  filter  can  not  adequately 
approximate  the  sharpness  of  the  true  peak  amplitude. 

The  random  number  generator  is  reset  a  total  of 
four  times,  producing  four  samples  at  each  noise  level 
for  each  depth.  After  the  statistical  data  are  obtained, 
a  linear  regression  analysis  of  the  four  data  sets  is 
made,  and  the  results  are  shown  in  Figure  16,  a  semi- 
logarithmic  plot  of  a  as  a  function  of  crossover 
wavenumber.  As  is  shown,  the  log  of  a  can  be 
represented  as  a  linear  function  of  kc  for  a  given  depth. 
The  slope  and  a-intercept  value  for  the  log10  of  a  (aQ) 
of  each  line  on  Figure  16  are  calculated  and  plotted  as  a 
function  of  depth  (Figures  17  and  18).  The  slopes  as  a 
function  of  depth  are  linear  as  are  log  of  aQ  as  a 
function  of  log  of  depth.  A  linear  regression  analysis 
of  these  two  factors,  slope  vs.  depth  and  log  of  aQ  vs. 
log  of  depth,  is  performed.  The  coefficients  of  the 
linear  regression  are  then  used  as  the  empirical 
relationship  for  a  and  kc,  given  as: 
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iogi0(a)  =  lQgio(ao)  +  s>kc 


S  =  -0.623  -  0.583*2 


log10(a0)  =  0.303  +  1. 928 • log10 (z) 


From  relationships  (a)  ,  (b)  ,  and  (c)  ,  an  optimal  value 


for  a  could  be  calculated  for  any  given  depth  and 


crossover  wavenumber. 


Since  only  a  small  range  of  depths  are  used  to 


establish  this  relationship,  some  error  might  be  expected 


for  larger  or  smaller  depths.  At  the  smaller  depths 


(0.0625  km  and  less),  the  linear  relationship  does  not 


strictly  hold,  but  this  may  be  due  to  a  lack  of  computer 


precision  with  a  values  on  the  order  of  10  and  less, 


Depths  greater  than  about  20.0  km  cannot  be  calculated 


due  to  machine  overflows  from  attempted  calculation  of 


the  exponential,  downward  continuation  function  at 


greater  depths. 


In  spite  of  a  lack  of  greater  depth  range,  it 


appears  that  a  is  not  a  highly  sensitive  parameter.  For 


example,  at  1.250  km,  standard  deviation  of  a  values  for 


a  given  noise  level  could  range  from  10  percent  of  the 


average  a  value  at  a  low  noise  level  to  50  percent  of  the 


average  a  value  at  a  higher  noise  level.  Despite  this 


high  variance,  rms  error  between  surface  data  regularized 


with  optimal  a  and  data  at  depth  dc  ranges  from  between  4 


percent  and  6  percent  of  the  maximum  value  of  the  true 


:  "V 


gravity  anomaly  at  depth,  the  error  level  increasing  with 
increasing  noise  level. 

Regularization  with  Optimal  Selection  of  a: 

The  empirical  relationship  between  a  and  kc 
described  in  the  last  section  is  tested  on  a  synthetic 
gravity  model  distinct  from  the  model  used  to  generate 
the  experimental  data.  The  gravity  effect  of  rectangular 
prismatic  sources  at  different  depths  is  calculated  at 
the  observation  plane,  z  =  0,  and  at  the  continuation 
depth,  dc,  and  the  procedure  of  Figure  19,  as  described 
below,  is  followed  to  test  the  accuracy  of  the  empirical 
a-kc  relationship.  For  each  depth,  width  of  the  source 
ranges  from  0.1  to  1.6  times  the  depth  of  the  source. 
Depth  extent  of  the  source  is  one  fourth  of  the  source 
depth  for  each  test  case.  Source  width  is  varied  in 
order  to  monitor  the  effect  of  finite  source  thickness  on 
efficiency  of  the  regularization  filter.  The  a  parameter 
produced  by  the  empirical  relationship  formula  compares 
favorably  with  the  "optimal"  a  value  calculated  using  the 
procedure  FMIN  as  previously  described.  Variations  on 
the  order  of  10  to  20  percent  are  observed  between  the  a 
values  using  the  two  different  methods. 

Rms  error  between  regularized,  downward  continued 
data  and  data  calculated  at  the  continuation  depth  ranges 
from  5  to  8  percent.  Again,  the  greatest  error  is 
observed  between  filtered  and  "true"  data  at  higher  noise 


*For  each  j,  the  depth  of  the  source  is  doubled;  the  range  of  depths  is  from 
0.125  to  18  km* 

For  k  =  1  to  m  do 

‘For  each  k,  the  width  of  the  source  is  doubled;  the  range  of  widths  is 
from  0.1  tc  1.6  times  the  depth  of  the  source.  The  depth  extent  of  the 
source  is  constant  at  0.25  times  the  depth  of  the  source* 

Generate  synthetic  gravity  data  at  the  surface  and  at  the  continuation 
depth  (0.9  times  the  modei  depth)  for  a  rectangular  source  at  depth  = 

2J ' 

For  i  =  1  to  4  do 

noise/signal  ratio  (n/s): 
if  i  =  1  then  n/s  =  2  percent 

if  i  =  2  then  n/s  =  5  percent 

if  i  =  3  then  n/s  =  10  percent 

if  i  =  4  then  n/s  =  20  percent 

*noise  to  signal  ratio  =  ratio  of  the  standard  deviation  of  the 
Gaussian  random  noise  to  the  maximum  gravity  anomaly  value* 

Add  noise  to  the  surface  data. 

Calculate  power  spectral  estimate  of  surface  data 

Find  kc  by  least  squares  spectral  mcaeling. 

Find  a  using  empirical  a  -kc  relationship  for  the  continuation 
depth  0.S*3j. 

Filter  the  data  using  the  regularization  ana  dcwnwara  continuation 
functions. 

Calculate  least  squares  error  between  the  filtered  surface  aata  ana 
data  at  the  continuation  depth. 

Compare  these  results  with  results  obtained  using  the  a  va.ue 
ootained  from  tne  ncminear  optimization  prcceaure 

Plot  results. 

End  for 

End  for 

Ena  for 

r.d 


Figure  19.  Procedure  for  testing  rectangular  source 
models  using  previously  derived,  empirical  a  vs.  kc 
relationship  as  a  function  of  depth. 
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levels.  Increasing  the  width  of  the  source  does  not 
increase  the  error;  however,  error  level  for  the 
rectangular  sources  is  higher  than  that  observed  for  line 
source  anomalies.  A  representative  example  of  one  of  the 
sources  used  (Figure  20)  and  results  for  two  different 
noise/signal  ratios  (0.05  and  0.20)  (Figures  21  and  22) 
are  presented. 

Testing  the  Forward  Solution: 

The  one-dimensional  Fourier  forward  gravity 
solution  with  the  addition  of  a  vertically  varying 
density  function  is  tested  on  a  synthetic  trapezoidal 
model  (Figure  23) ;  the  density  function  is  arbitrarily 
chosen  and  specified  a  priori .  Results  of  the  forward 
solution  are  compared  to  results  obtained  from  a  constant 
density,  forward  gravity  technique  which  sums  the  gravity 
effects  of  prismatic  bodies  (Corbato,  1965) .  A  linear 
density  function  is  simulated  for  the  constant  density 
solution  by  calculating  the  gravity  effect  of  multiple 
layers.  Density  contrast  at  each  layer  is  equal  to  the 
value  of  the  specified  density  function  at  the  midpoint 
of  each  layer.  This  corresponds  to  a  straight  numerical 
solution  of  the  field  integral.  Increasing  the  number  of 
layers  alters  the  results  only  slightly  and  does  not 
lessen  variation  of  the  gravity  anomaly  values  between 
the  two  methods;  therefore,  resolution  of  the  field 
integral  method  is  not  increased  as  a  result  of 
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Figure  20.  Example  of  rectangular  source  model  at  6  km 
depth  and  gravity  effects  calculated  at  the  surface  or 
observation  plane  (z  =  0)  and  at  the  continuation  depth 


(z  =  dr)  where  dr  =  0.9  times  the  source  depth. 
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Figure  21.  Example  of  results  of  applying  the 
procedure  described  in  Figure  19  for  the  source  in 
Figure  20  using  empirical  a  vs.  kc  relationship  as  a 
function  of  depth:  (a)  surface  data  with  2  percent 
noise/signal  ratio,  (b)  spectral  model  with  kc  at  0.679 
radians/km,  (c)  regularization  filter  plotted  with  the 
data  spectrum  and  a  =  0.143,  (d)  comparison  of 
regularized,  downward  continued  surface  data  with  data 
calculated  at  the  continuation  depth  (5.4  km)  with  an 
rms  error  of  0.8016  mGals. 
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Figure  22.  Results  using  the  same  model  as  in  Figure 
21  with  noise/signal  ratio  of  20  percent:  (a)  surface 
data  with  added  noise,  (b)  spectral  model  with  kc  at 
0.405  radians/km,  (c)  regularization  filter  plotted 
with  data  spectrum  and  a  =  1.547,  (d)  comparison  of 
regularized,  downward  continued  surface  data  with  data 
calculated  at  the  continuation  depth  (5.4  km)  with  an 
rms  error  of  0.8965  mGals. 
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increasing  the  parameterization.  Results  presented  here 


are  calculated  with  a  12  layer  model.  The  Fourier 


solution  is  compared  with  summed  results  of  the  field 


integral  solution.  Rms  error  between  the  two  results  is 


1.43  mGals  (4.1  percent  error)  and  maximum  error  is  3.25 


mGals  (9.4  percent  error)  with  greatest  error  occurring 


along  the  flanks  of  the  anomaly  (Figure  24) .  Discrepancy 


between  the  two  methods  can  be  attributed  to  a  difference 


in  source  parameterization.  For  the  field  integral 


method,  each  layer  is  parameterized  with  four  vertices 


and  a  single  density  contrast.  Model  parameterization 


for  the  Fourier  calculation  assumes  a  continuous  body 


sampled  at  discrete  intervals  with  a  smoothly  varying, 


linear  density  function. 


The  difference  is  most 


pronounced  in  its  affect  at  the  shallow  body  corners. 


Testing  the  Inverse  Solution: 


The  inverse  solution  is  tested  using  a  synthetic, 


basin-type  gravity  model.  Gaussian  random  noise  is  added 


to  gravity  anomaly  values  calculated  from  the  synthetic 


model.  In  addition,  data  are  pre-processed  by  applying  a 


5  percent  cosine  taper  and  padding  the  data  with  zeros 


out  to  twice  the  interval  length.  An  initial  guess  for 


the  model  is  made  using  a  Bcuguer  slab  approximation, 


di  =  Sg!  ,  (2t  G  c] 


where : 


=  slab  thickness 
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Figure  24.  Results  of  testing  the  Fourier  forward 


gravity  algorithm  against  the  field  integral  solution 


of  Corbato  (1965)  for  the  model  of  Figure  23 


differences  along  the  sides  of  the  anomaly  are  probably 


due  to  differences  in  model  parameterization 
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a  =  density  contrast 

Sg^  =  pre-processed  gravity  data 

and, 

i  =  1  to  N,  where  N  =  length  of  data  interval. 

A  first  approximation  to  the  continuation  depth 
used  in  the  choice  for  a  and  in  the  regularization  filter 
is  a  depth  halfway  between  minimum  and  maximum  model 
depths  estimated  by  the  Bouguer  slab  formula  for  the 
initial  guess.  Density  contrast  used  in  the  initial 
guess  is  specified  a  priori  and  is  the  value  of  the 
specified  density  function  midway  between  a  minimum  and 
maximum  model  depth  estimate. 

The  inverse  modeling  procedure  as  adapted  from 
Ferguson  et  al .  (1988)  is  shown  in  Figure  25.  The  data 

are  regularized  using  the  empirical  a-kc  relationship, 
with  crossover  wavenumber  determined  by  the  spectral 
modeling  procedure  previously  discussed  and  continuation 
depth  estimated  from  the  initial  guess.  After  each 
iteration  of  the  inversion  algorithm,  a  new  model  is 
generated  and  used  in  the  forward  algorithm  to  generate 
synthetic  gravity  data.  Synthetic  data  are  compared  with 
observational  data  and  rate  of  change  of  the  ratio  of  rms 
error  between  the  two  data  sets  from  the  previous  to 
present  solution  is  used  to  determine  convergence  of  the 
solution.  If  convergence  is  not  obtained,  the 

continuation  depth  of  the  present  model  is  used  to 


Begin 

Read  data  and  compute  initial  guess  model. 

Take  Fourier  transform  of  data. 

Calculate  power  spectral  estimate  of  data. 

Find  crossover  frequency  (kc)  using  least  squares  spectral  modeling. 

Find  regularization  parameter,  a  .using  empirical  a  vs.  kc  relationship  as  a 
function  of  continuation  depth  (zc). 

Perform  one  evaluation  of  Fourier  forward  gravity  algorithm  using  initial  guess 
model  and  calculate  error  between  theoretical  and  observed  data  to  set  initial  rms 
error  level. 

While  the  number  of  iterations  is  less  than  the  maximum  (specified  a  priori  1  and 
the  ratio  of  current  to  previous  least  squares  error  is  greater  or  equal  to  the 
tolerance  level  (specified  a  priori)  do 

Take  Fourier  transform  of  data. 

Filter  the  data  using  the  regularization  function. 

Invert  the  data  using  the  Fourier  gravity  inversion  aigoritnm. 

Find  new  a  using  empirical  relationship  at  new  zc. 

Filter  the  model  using  the  regularization  function. 

Calculate  theoretical  gravity  data  by  using  the  new  mcaei  in  an  aeration  of 
the  Fourier  forward  gravity  algorithm. 

Find  the  least  squares  error  between  the  observed  ana  theoretical  aala 
Calculate  ratio  of  current  to  previous  least  squares  error 
End  while 
Plot  the  results. 


Figure  25.  Procedure  for  inversion  of  gravity  data 
using  regularization  and  the  Fourier  gravity  inversion 
technique  with  a  linear,  vertically  varying  density 


function . 


regularize  the  original  data  and  another  iteration  is 
then  performed. 

In  practice,  rms  error  is  "windowed"  so  as  not  to 
extend  beyond  the  length  of  the  original  interval  of 
observations.  Windowing  the  error  speeds  convergence  due 
to  exclusion  of  error  from  any  Gibbs  oscillation  along 
the  edges  of  the  anomaly,  that  error  having  greater 
weight  than  error  due  to  a  misfit  of  data  representing 
deeper  portions  of  the  model.  Differential  error  weights 
result  since  density  is  a  function  of  depth.  Resolution 
at  greater  depths  is  less  as  the  value  of  the  density 
contrast  is  less  (for  a  negative  density  contrast) 
resulting  in  a  disparity  in  the  contribution  to  the 
overall  gravity  effect  between  deep  and  shallow 
components.  Without  windowing,  this  disparity  tends  to 
bias  the  inversion  procedure  toward  a  better  fit  of  the 
data  on  the  sides  of  the  anomaly  rather  than  on  the  peak. 

Convergence  criteria  are  determined  a  priori .  In 
general,  a  maximum  of  ten  iterations  are  sufficient  and 
rate  of  change  of  the  ratio  of  least  squares  error 
between  synthetic  and  calculated  gravity  anomaly  values 
of  between  one  and  five  percent  is  considered 
satisfactory.  Speed  of  convergence  depends  on  a 
reasonable  initial  guess  of  the  gravity  model;  a  good 
initial  guess  reduces  the  number  of  iterations  required 
for  convergence.  Underestimating  the  density  contrast 


used  in  initial  depth  estimates  tends  to  produce  models 
which  underestimate  depth  of  the  source  and  converges 
more  slowly.  Overestimation  of  the  density  contrast 
results  .in  more  rapid  convergence;  however,  if  model 
depths  obtained  from  initial  guess  estimates  occur  at 
depths  which  exceed  the  saturation  level  of  the  density 
function,  the  solution  diverges. 

For  synthetic  model  results,  an  rms  error  of  2.52 
mGals  or  about  3  percent  of  the  maximum  gravity  anomaly 
value  is  observed;  however,  effects  of  the  addition  of 
noise  with  a  2  percent  noise/signal  ratio  are  completely 
filtered  out  (Figure  26)  .  Rms  error  between  the  gravity 
model  obtained  from  inversion  and  the  synthetic  model 
(Figure  27)  is  0.087  km  or  about  2  percent  of  the  maximum 
depth  whereas  maximum  error  is  0.18  km  or  about  4.5 
percent  error. 

For  comparison,  the  same  synthetic  model  is  used 
and  a  noise/signal  ratio  5  percent  is  added  to  the 
synthetic  gravity  data.  The  final  anomaly  fits  the 
synthetic  anomaly  very  well  in  a  least  squares  sense, 
with  an  rms  error  of  4.3  mGals  or  5  percent  of  the 
maximum  anomaly  value  (Figure  28)  .  Model  depths  also 
compare  favorably  at  the  higher  noise  level  (Figure  29) 
with  an  rms  error  of  0.11  km  or  2.7  percent  and  maximum 
error  of  0.25  km  or  6.3  percent.  Greatest  error  in  both 
cases  occurs  along  the  sides  of  the  model,  but  error  in 


Figure  26.  Comparison  of  theoretical  data  calculated 


from  gravity  model  obtained  from  the  inversion 
procedure  of  Figure  25  compared  with  gravity  data 
calculated  froma  synthetic  basin  type  model  with  2 
percent  noise/signal  ratio. 
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Figure  28.  Comparison  of  theoretical  data  calculated 


from  gravity  model  obtained  from  the  inversion 


procedure  of  Figure  25  compared  with  gravity  data 


calculated  from  a  synthetic  basin  type  model  with  5 


percent  noise/signal  ratio. 
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both  cases  is  considered  well  within  acceptable  error 
bounds . 

Error  level  can  be  decomposed  into  contributions 
from  observational  error,  the  added  noise  level,  and 
resolving  error,  resulting  from  suppression  of  high 
wavenumber  information  during  regularization.  Since  the 
regularization  filter  is  optimized  based  on  the 
assumption  of  sources  at  depth,  resolution  of  the  model 
along  the  shallow  flanks  of  is  poorer  as  the 
corresponding  short  wavelength  information  is  filtered 
out.  Width  of  the  regularization  filter  is  a  qualitative 
measure  of  model  resolution  (Ferguson  et  al.,  1988). 
Since  for  both  noise  levels  crossover  wavenumber  is 
relatively  small  (Figure  30) ,  bandwidth  of  the  low  pass, 
regularization  filter  is  also  small  resulting  in  a  lack 
of  resolution  as  observed  along  the  edges  of  the 
synthetic  basin-type  model. 


Collection  and  Reduction  of  Gravity  Data 
at  Pahute  Mesa .  Nevada 

The  data  used  in  this  study  consist  of 
approximately  1700  gravity  stations  collected  by  the  USGS 
and  the  University  of  Texas  at  Dallas  (UTD) . 
Approximately  250  gravity  stations  and  elevations  were 
collected  by  UTD  in  June,  1985  at  Pahute  Mesa,  Nevada 
along  detailed  survey  lines  (station  spacing  of 
approximately  0.3  km).  The  gravity  data  were  collected 
with  two  LaCoste  and  Romberg  gravimeters  and  elevations 
were  shot  in  at  each  station  with  a  Zeiss  Elta-4 
electronic  distance  meter.  The  accuracy  of  these 
stations  is  probably  good  to  within  +0.5  raGals.  The 
detailed  gravity  stations  were  merged  with  the  existing 
approximately  1400  USGS  gravity  stations  at  Pahute  Mesa 
(Figure  31) .  The  accuracy  of  the  existing  data  is 
approximately  +5  raGals  (Healy,  1968) .  After  merging  the 
two  data  sets,  the  data  were  edited  and  stations  deleted 
that  were  considered  bad  points. 

The  complete  Bouguer  anomaly  map  of  Pahute  Mesa 
(Figure  32)  was  obtained  by  the  following  reduction 
procedure.  Hand  terrain  corrections  of  Hammer  zones  B-D 
(Hammer,  1939)  were  calculated  for  both  data  sets  and 
added  to  outer  zone  terrain  corrections  calculated  from 
the  terrain  correction  program  of  Plouff  (1977).  Total 
terrain  corrections  ranged  from  about  0.2  to  2.0  mGals. 
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Reduction  density  for  the  complete  Bouguer  anomaly  was 
2000  kg/m  .  This  density  was  an  average  density  of  the 
rocks  constituting  the  topography  at  Pahute  Mesa  down  to 
a  depth  of  about  0.15  km,  which  ranged  from  about  1600  to 
2200  kg/m  ,  estimated  from  analysis  of  average  densities 
obtained  from  borehole  gravity  logs.  The  data  were 
gridded  at  0.5  km  intervals  and  contoured  at  1  mGal. 

Possible  error  in  the  UTD  data  could  be  attributed 
to  a  number  of  factors.  A  list  of  some  of  these  factors 
is  given  below. 

1.  Equipment  failure .  Equipment  failed  on  several 
occasions  due  to  afternoon  high  temperatures  at  Pahute 
Mesa,  sometimes  well  over  100  F.  Parts  of  the  survey 
where  this  occurred  were  redone;  however,  the  possibility 
of  error  remains. 

2.  Observational  error.  Repeatability  to  within 
0.01  mGals  was  required  for  each  gravity  station,  to 
within  ±1.5  meters  for  vertical  distances,  and  to  within 
±0.030  kilometers  for  horizontal  distances.  Some  error 
is  to  be  expected  from  inaccurate  posting  of  the  data  or 
inaccurate  reading  of  instruments. 

3.  Plotting  error.  When  converting  digital  survey 
data  to  map  data,  some  plotting  errors  can  be  expected. 
In  this  case,  the  data  was  plotted  by  hand  onto  a  base 
map  at  a  scale  of  1:24000,  and  the  x-y  locations  for  each 
station  wnre  digitized  from  the  base  map. 


Horizontal 


locations  at  this  scale  are  only  good  to  about  +0.030 
kilometers . 
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A  Pahute  Mesa  Gravity  Model 
The  inversion  procedure  is  tested  next  on 
observational  gravity  data  from  Pahute  Mesa,  Nevada.  A 
complete  Bouguer  anomaly  map  of  the  area,  contoured  at 


one  mGal  and  with  a  reduction  density  of  2000  kg/m  ,  is 
shown  in  Figure  32  with  model  profiles  A-A'  and  B-B' 
indicated.  Station  locations  are  shown  in  Figure  31. 
Only  profile  B-B'  will  be  discussed  in  this  report,  due 
to  unresolved  stratigraphic  details  along  A-A' . 

Profile  locations  are  chosen  parallel  (A-A')  and 
perpendicular  (B-B')  to  the  northeast-trending, 
elliptically-shaped  gravity  low  over  Pahute  Mesa  thought 
to  outline  Silent  Canyon  caldera  (Healy,  1968;  Orkild  et 
al.,  1968).  The  caldera  rim,  as  outlined  by  steepening 
of  the  gravity  contours,  is  obscured  to  the  southwest  due 
to  thickening  of  Timber  Mountain  volcanics  (Orkild  et 
al.,  1968).  Based  on  the  density  structure  of  Pahute 
Mesa,  Silent  Canyon  caldera  can  be  modeled  as  a  higher- 
density  half  space  corresponding  to  the  Belted  Range  Tuff 
below  a  lower  density  layer,  the  caldera  fill. 
Therefore,  this  is  essentially  a  basin-type  model.  Based 
on  this  assumption,  choice  of  the  Fourier  gravity 
inversion  technique  seems  appropriate  as  the  model 
parameters  produced  by  the  inversion  will  delineate  a 
detailed  model  of  the  caldera  at  the  interface  between 


the  Belted  Range  tuff  and  the  caldera  fill.  The  addition 
of  a  vertically  varying,  linear  density  function  makes  it 
possible  to  fit  the  linear  density  vs.  depth  relationship 
as  discussed  earlier  into  the  model. 

In  order  to  constrain  the  regional  anomaly, 
profiles  have  been  extended  outside  the  edges  of  the 
caldera  to  gravity  stations  resting  on  outcrops  of  rocks 
representing  the  gravity  basement,  i.e.,  the  Belted  Range 
tuff  or  an  equivalent  lithology.  Profile  B-B'  is 
constructed  to  extend  northwest  from  outcrop  of  the 
Grouse  Canyon  member  of  the  Belted  Range  Tuff  southeast 
to  outcrop  of  Quartet  Dome  lavas.  Quartet  Dome  lavas 
are  lithologically  equivalent  to  (i.e.,  have  about  the 
same  density  contrast  as)  and  are  found  interbedded  with 
Grouse  Canyon  tuffs. 

The  observed  anomaly  is  assumed  to  be  the  sum  of 
two  components,  (1)  a  regional  anomaly  of  wavelength 
equal  to  or  greater  than  the  profile  length  and  (2)  a 
residual  anomaly  composed  of  smaller  wavelength  effects 
(Nettleton,  1976)  .  Based  on  the  assumption  that  the 
Pahute  Mesa  gravity  low  is  caused  by  lower  density 
caldera  fill,  the  residual  gravity  anomaly  at  basement 
outcrops  at  either  end  of  the  profiles  should  be  zero. 
Therefore,  the  regional  anomaly  at  either  profile 
terminus  is  defined  by  the  observed  Bouguer  anomaly  value 
at  that  location.  Shape  of  the  regional  anomaly  is 


assumed  to  be  a  low  order  surface  with  wavelength 
proportional  to  the  profile  length.  The  regional  anomaly 
for  Pahute  Mesa  is  fitted  with  a  second  degree 
polynomial . 

The  regional  anomaly  calculation  is  incorporated 
into  the  modeling  procedure  using  a  method  adapted  from 
Ferguson  et  al.  (1988)  (Figure  33).  Boreholes  in  Pahute 
Mesa  which  penetrate  the  basement  along  each  profile  line 
are  used  to  constrain  the  regional  (Figure  34)  .  As 
outlined  in  Figure  33,  after  each  gravity  inversion,  a 
new  model  is  produced.  Model  depths  are  compared  to 
borehole  depths  and  residual  errors  used  to  calculate  a 
new  regional  anomaly.  After  removal  of  the  new  regional, 
another  new  model  is  calculated.  The  modeling  process  is 
terminated  when  residual  error  between  borehole  depths  to 
the  basement  interface  and  model  depths  is  within 
acceptable  limits. 

In  addition  to  removing  the  regional  anomaly,  gravity 
effects  of  the  higher  density  Lower  Rainier  Mesa  tuffs 
and  Area  20  lavas  are  stripped  from  the  residual  anomaly 
based  on  thickness  estimates  from  borehole  stratigraphic 
data.  Average  density  contrasts  of  400  kg/m3  and 
250  kg/m3  are  assumed  for  Rainier  Mesa  tuffs  and  Area  20 
lavas ,  respectively,  based  on  analysis  of  the  density  vs . 
depth  plots  from  the  borehole  gravity  logs  (Figure  5)  . 
The  gravity  effects  of  these  layers  are  removed  using  a 


Begin 

Fit  smoothing  spline  to  profile  gravity  data  and  interpolate  at  desired  sample 
spacing. 

Fit  smoothing  spiine  to  gravity  effects  of  higher  density  lavas  and  tuffs 
calculated  along  model  profiles  and  subtract  from  profile  observations. 

Find  boreholes  along  profile  line  which  penetrate  to  the  gravity  basement  to  use 
as  constraints  for  the  regional  anomaly  calculation. 

Extend  model  profile  to  include  gravity  stations  resting  on  outcrop  of  lithology 
representing  gravity  basement  for  constraint  of  regional  anomaly. 

Estimate  initial  regional  anomaly  by  a  straight  line  fit  of  gravity  stations  on 
basement  outcrops. 

Subtract  initial  regional  anomaly  from  observed  gravity  anomaly. 

Execute  inversion  algorithm  (Figure  25)  to  produce  an  initial  gravity  model. 
Calculate  residual  error  between  borehole  depths  and  initial  model  depths. 

Repeat 

Use  residual  error  between  borehole  depths  and  model  depths  to  update 
gravity  values  used  in  the  regional  anomaly  calculation. 

Ft  a  second  degree  polynomial  by  least  squares  to  updated  gravity  values  for 
new  regional  anomaly  estimate. 

Subtract  regional  anomaly  from  observed  gravity  anomaly. 

Execute  inversion  algorithm  (Figure  25)  to  produce  a  new  gravity  model. 

Calculate  residual  error  between  borehole  depths  and  modei  depths. 

Until  residual  error  between  borehole  depths  and  model  depths  is  less  than 
desired  tolerance. 

Interpret  mcdel. 


Figure  33.  Algorithm  outlining  the  modeling  procedure. 
Several  different  programs  are  represented  by  this 
algorithm. 
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Bouguer  slab  approximation.  The  residual  data  are  fitted 
with  a  smoothing  spline  (Reinsch,  1967)  and  interpolated 
at  0.5  km  intervals.  Approximately  57  gravity  stations 
are  used  for  profile  B-B'. 

Figure  35  is  the  gravity  model  for  profile  B-B’  . 
This  profile  is  of  particular  importance  as  it  closely 
parallels  the  profile  of  a  seismic  reflection/refraction 
survey  which  was  run  in  June,  1986  at  Pahute  Mesa.  Least 
squares  error  of  the  theoretical  gravity  data  is 
+1.27  mGals  (4  percent  error)  with  a  maximum  error  of 
6.78  mGals  (21  percent  error)  (Figure  36).  Rms  error  of 
the  depth  residuals  is  +0.048  km  with  a  maximum  of  0.16 
km  at  borehole  Ue20e-1.  The  regional  anomaly  is  shown  in 
Figure  37.  The  spectral  model  for  profile  B-B'  is  shown 
in  Figure  38. 

Effectiveness  of  the  Regularization  Filter: 

As  the  results  have  demonstrated,  application  of  a 
regularization  filter  stabilizes  downward  continuation  by 
suppressing  high  wavenumber  components  of  the  data 
spectrum.  A  basic  assumption  of  the  experimental  results 
is  that  the  gravity  anomaly  was  composed  of  sources  at 
depth  and  that  effects  of  any  intervening  masses  are 
filtered  out  in  order  to  stabilize  the  downward 
continuation  process.  If  the  gravity  anomaly  is  due  to 
contributions  from  a  source  at  a  single  interface  or 
depth,  this  assumption  is  valid.  However,  since 
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Figure  36.  Observed  residual  gravity  anomaly  plotted 
with  gravity  anomaly  calculated  in  the  modeling 
procedure  (Figure  33)  for  profile  B-B’ . 
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Figure  37.  Observed  Bouguer  anomaly  values  for  profile 


B-B'  plotted  with  the  final 
regional  anomaly  for  profile 
polynomial . 


regional  anomaly.  The 
B-3 1  is  a  second  degree 


Figure  38.  Spectral  model  for  profile  B-B' .  Profile 
B-B'  appears  to  have  shallow  sources  which  have  been 
filtered  out  in  addition  to  the  noise. 
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contributions  to  the  gravity  effect  from  shallower 
sources  are  suppressed  in  the  filtering  process,  the 
resulting  gravity  model  from  data  containing 
contributions  from  sources  at  different  depths  will  not 
contain  any  information  about  shallower  sources.  For 
example,  for  Pahute  Mesa  gravity  models,  greatest  data 
misfit  for  each  profile  occurs  along  the  sides  of  the 
anomaly  where  the  caldera  depth  shallows.  This  is  due  to 
suppression  of  high  wavenumber  content  of  the  gravity 
data  during  the  regularization  process.  An  improvement 
to  the  regularization  process  would  be  to  allow  for  non¬ 
stationary  aspects  of  gravity  data  for  areas  such  as  a 
basin,  where  basement  depth  varies  as  a  function  of 
distance  across  the  basin.  This  could  be  accomplished  by 
assuming  that  the  data  spectrum  varies  smoothly  over 
space,  and  spectral  estimates  could  be  made  locally  to 
accommodate  changes  in  depth  across  an  area  (Priestly, 
1981) .  By  using  a  space-variable  regularization  filter, 
effects  of  shallow  structure  could  be  preserved  and 
filter  width  used  as  an  estimate  of  resolution. 
Covariance  of  the  model  could  also  be  estimated  from  the 
filter  function. 

Effectiveness  of  the  Vertically  Varying  Density  Function: 

Although  effective  at  Pahute  Mesa,  use  of  a 
linearly  varying  density-depth  relation  has  serious 
drawbacks  which  limiting  the  overall  usefulness  of  such  a 


method.  In  the  procedure  described  in  Figure  33,  the 
inverse  iteration  is  repeated  until  convergence  is 
obtained,  as  previously  discussed.  If,  during  execution 
of  the  procedure,  model  depth  falls  below  the  depth  at 
which  the  value  of  the  density  function  equals  zero,  the 
solution  begins  to  diverge  rapidly.  Model  depths  could 
be  constrained  to  have  a  lower  bound  set  above  turnover 
point  of  the  density  function.  However,  this  is  not  a 
totally  satisfactory  solution  as  it  biases  the  model  with 
respect  to  the  imposed  constraint  (Menke,  1984) .  For 
Pahute  Mesa,  the  density  function  saturates  at  a  depth  of 
about  2.0  km,  and  parts  of  the  caldera  basement  which  are 
deeper  than  2.0  km  cannot  be  modeled  using  the  density 
function  presented  in  this  study.  Since  the  density 
function  adopted  is  somewhat  arbitary,  based  on  a  linear 
fit  to  boreholes  with  density  information,  changing  the 
density  function  to  accommodate  deeper  parts  of  the 
caldera  will  circumvent  saturation  of  the  density 
function.  However,  a  better  solution  is  to  use  the  more 
appropriate  exponential  density  function  which  allows  for 
falloff  rate  and  saturation  of  density  at  depth.  The 
linear  density  variation  of  the  Fourier  gravity  inversion 
method  is  then  most  useful  when  the  density  function 
defining  the  model  does  not  change  sign  near  the  basement 


interface . 


Estimates  of  covariance  and  resolution  for  non¬ 
linear  inverse  problems  are  net  easily  calculated  since 
the  distribution  function  is  not  usually  known  (Menke, 
1984)  .  In  the  case  of  the  Fourier  gravity  inversion 
method  with  the  addition  of  a  linearly  varying  density 
vs.  depth  function,  it  is  observed  that  resolution 
decreases  with  increasing  depth  as  density  contrast 
becomes  smaller.  Perturbations  to  the  model  at  shallower 
depths  contributes  more  significantly  to  the  overall 
gravity  effect  than  perturbations  at  greater  depths. 
Weighting  least  squares  error  calculations  as  a  function 
of  depth  could  compensate  for  this  lack  of  resolution. 
Also,  as  the  regularization  filter  suppresses  effects  of 
shallower  sources,  model  resolution  at  greater  depths  is 
decreased  in  order  to  decrease  model  variance. 

Another  problem  with  the  Fourier  gravity  inversion 
technique  is  the  tendency  for  model  depths  along  profile 
edges  to  become  negative,  which  is  a  result  of 
overfiltering  along  the  flanks  of  the  caldera,  as 
discussed  previously.  Because  of  the  steep  gravity 
gradient  along  profile  edges,  the  loss  of  information 
from  regularization  results  in  a  lack  of  continuity  along 
the  profile  edges.  According  to  Dugue's  Theorem,  the 
Fourier  transform  of  a  continuous,  positive  definite 
function  is  uniformly  convergent;  however,  since  the 
discrete  representation  of  positive  definite  functions 


are  generally  not  continuous  at  the  endpoints,  uniform 
convergence  is  not  assured  (Schmetterer ,  1965) . 
Interpretation  of  Pahute  Mesa  Depth  Models: 

Choice  of  the  inverse  Fourier  modeling  technique 
for  Pahute  Mesa  is  made  based  on  the  assumption  that 
fluctuations  in  the  gravity  anomaly  are  caused  by 
variations  in  the  source  topography.  Suppression  of 
shallower  source  effects  is  desirable  to  decrease  model 
variance.  However,  since  the  assumption  is  also  made  to 
model  Silent  Canyon  Caldera  as  a  basin-type  structure, 
inclusion  of  a  space-variable,  regularization  filter 
would  increase  model  resolution  along  the  edges  of  the 
caldera  where  signficant  structural  changes  are  apparent 
and  the  overfiltering  in  this  study  occurs. 

For  example,  in  the  gravity  model  for  profile  B- 
B',  the  most  significant  mismatch  between  observed  and 
calculated  anomalies  is  toward  the  southeastern  profile 
edge  at  a  depth  of  about  0.9  km.  This  shallow  structure 
coincides  with  a  mapped  surface  fault,  the  Scrugham  Peak 
fault,  (Orkild  et  al.,  1969)  which  appears  from  this 
model  to  extend  into  the  basement.  The  spectral  model 
for  profile  B-B'  indicates  the  presence  of  at  least  one 
and  possibly  two  shallow  sources.  The  poor  least  squares 
fit  to  the  periodogram  results  in  larger  maximum  error 
between  the  residual  anomalies  as  high  wavenumber  content 
of  the  anomaly  corresponding  to  shallower  structure  is 
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filtered  out.  Model  resolution  and  variance  for  this 
profile  can  be  improved  by  use  of  a  space  variable, 
regularization  filter. 

Profile  B-B'  is  well  constrained  by  boreholes 
(Figure  34) ,  indicating  that  basement  depths  are  probably 
good  to  within  the  least  sguares  residual  error  between 
model  and  borehole  depths.  Also,  the  good  distribution 
of  boreholes  provides  better  control  on  thicknesses  of 
higher  density  units  stripped  from  profile  data.  Because 
of  the  abundance  of  normal  faulting  at  Pahute  Mesa, 
thickness  of  the  higher  density  units  varies  considerably 
in  some  places. 

Profile  B-B'  does  not  transect  the  center  of 
gravity  low;  however,  undulations  in  the  deepest  part  of 
the  model  for  B-B'  suggest  some  structural  complexity  in 
basement  topography  of  the  caldera.  A  three-dimensional 
model  of  Pahute  Mesa  data  could  resolve  some  of  these 
questions . 

The  shape  of  the  second  degree  polynomial  regional 
anomaly  subtracted  from  each  profile,  a  regional  high 
indicating  a  positive  density  contrast,  may  indicate  an 
intrusive  source  at  depth.  The  existence  of  a  large 
intrusive  body  has  been  suggested  by  Snyder  and  Carr 
(1984)  in  a  study  area  southwest  of  Timber  Mountain  based 
on  aeromagnetic,  borehole,  heat  flow  and  gravity  data. 
Spence  (1974)  postulated  the  existence  of  a  high  velocity 


zone  in  the  upper  mantle  beneath  Silent  Canyon  caldera 
from  analysis  of  teleseismic  P-wave  residuals.  If  this 
zone  of  high  velocity  corresponds  to  a  crustal  or  upper 
mantle  magma  body  as  suggested  by  Spence,  then  the 
gravity  data  should  exhibit  an  associated  low  order,  long 
wavelength  gravity  high  unless  it  is  isostatically 
compensated. 


The  Wide  Angle  Seismic  Survey 


In  the  specification  of  a  geophysical  model  at 
Pahute  Mesa  the  borehole  data  are  limited  in  depth  extent 
and  coverage  and  the  gravity  models  lack  resolution  at 
depth  and  are  not  directly  sensitive  to  seismic  velocity. 
Direct  seismic  observation  of  the  caldera  structure  would 
be  of  considerable  importance  to  the  creation  of  a 
geophysical  model.  The  presence  of  local,  thin,  high 
velocity  layers  (such  as  the  Rainer  Mesa  tuff)  and  a  very 
low  Q  in  the  near  surface  materials  (perhaps  as  low  as 
30)  make  conventional  seismic  reflection  surveys 
impractical.  With  a  plan  based  on  some  preliminary 
hypotheses  concerning  the  caldera  structure,  an 
unconventional,  wide  angle  seismic  survey  has  been 
performed.  In  addition  to  the  geophysical  constraints  the 
presence  of  archeological  sites  on  Pahute  Mesa  limit  the 
possible  location  of  shot  holes. 

Geologic  and  gravity  models  suggest  that  the  Silent 
Canyon  caldera  can  be  thought  of  as  a  basin  founded  on 
older,  less  porous,  peralkaline  volcanics  associated  with 
the  caldera  collapse  (termed  "basement"  in  this  report) . 
The  basin  is  filled  with  low  density,  high  porosity  tuffs 
and  higher  density  lava  and  tuff  flows  from  sources  to 
the  south  and  west.  The  basement  rocks  outcrop  to  the 
southeast  of  Pahute  Mesa.  The  caldera  margin  could  be 
steeply  faulted  or  more  gently  slumped  with  large  zones 
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of  brecciation.  Several  major  north-south  trending  faults 
cut  through  the  area  and  have  been  active  during  much  of 
the  volcanic  deposition.  A  major  northwest  trending 
structure  may  also  be  present. 

A  joint  program  for  the  execution  of  a  seismic 
experiment  was  developed  between  UTD  and  Los  Alamos  in 
1985  and  1986.  Richard  Warren  and  Allen  Cogbill  of  the 
Earth  and  Space  Sciences  Division  at  Los  Alamos  were 
active  in  the  planning  of  the  experiment.  Los  Alamos 
provided  field  support  at  NTS,  field  personnel,  some 
equipment  and  shot  hole  drilling.  Ed  Criley  of  the  USGS 
at  Menlo  Park,  CA  was  subcontracted  to  acquire  and  handle 
the  explosives.  Based  on  the  considerations  discussed 
above  the  program  illustrated  in  Figure  39  was  laid  out 
in  the  Spring  of  1986.  A  single  shot  point  was  located  on 
outcropping  basement  volcanics  in  Big  Burn  Valley 
southeast  of  Pahute  Mesa.  A  line  of  recording  sites  was 
drawn  along  available  roads  on  the  Mesa  above  and 
trending  toward  the  hole  PM-1.  This  hole  is  one  of  the 
deepest  drilled  at  Pahute  Mesa  and  is  a  logical  place  to 
tie  the  end  of  the  line.  The  line  is  close  to  several 
other  boreholes  providing  stratigraphic  and  velocity 
control.  The  hole  U19w  is  of  particular  interest  as  it 
appears  to  be  near  the  caldera  margin.  Both  the  caldera 
margin  (as  indicated  on  Orkild  et  al.,  1969)  and  the 
north-south  trending  faults  are  crossed  at  a  favorable 


91 


angle  by  the  line.  The  profile  is  unreversed  by  shots  at 
the  northwestern  end  of  the  line  and  major  lateral 
structural  variations  occur  along  the  line.  The  intent  is 
to  inject  seismic  energy  directly  into  the  basement  and 
observe  refracted  and  reflected  waves  from  the  basement 
after  propagation  through  the  post  caldera  structure 
above.  The  interpretation  of  undershot  structures  similar 
to  this  was  successful  at  Crater  Flat  to  the  south  of 
this  area  (Hoffman  and  Mooney,  1984) .  The  interpretation 
will  rely  on  geologic  and  gravity  models  to  insure  some 
measure  of  uniqueness.  One  goal  of  this  experiment  is  to 
demonstrate  that  useful  seismic  data  can  be  acquired  at 
Pahute  Mesa.  If  this  turns  out  to  be  the  case,  future 
experiments  will  certainly  be  planned  with  more  complete 
shot  point  distributions. 

The  recording  line  is  12.5  km  long  with  a  3.5  km 
offset  from  the  shot  point,  necessitated  by  topography. 
Twenty  four  channel  spreads  of  1.5  km  aperture  and  67  m 
group  interval  were  deployed.  Each  channel  recorded  a 
string  of  six  4.5  Hz  geophones  spaced  3  to  6  m  apart.  A 
Geometries  ES-2420  recording  system  was  used  for  15  b:*: 
floating  point  digital  records  on  nine  track  tape.  A  . 
msec  sample  interval  was  used.  The  spreads  were  shot  *•  - 
the  Big  Burn  Valley  shot  point  and  also  :r'~  i  ■* 
with  a  12  liter  Dinose  is  ( ARCO  trademark  o  i  o 
The  spreads  recording  the  explosive  shots  r  i  • 


»v  r 


UNCLASSIFIED 
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Flag  #'s:  Sj  =  flag  locations  near  the  HE-shot  point  (Big  Burn  Valley) 
All  other  flag  #‘s  represent  locations  on  Pahute  Mesa 

SP  »  flag  #  of  shotpoint 

NEAR  OFFSET  *  offset  from  high  explosive  shotpoint  in  Big  Burn  Valley 
NVSS  =  #  of  vertically  sunned  shots  yielding  this  record 
FILL  TM  =  fill  time  of  Dinoseis  gins  in  seconds 
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Figure  40.  The  data  from  the  1986  seismic  survey  is 
shown  in  record  section  form.  The  data  have  been 
lowpass  filtered  at  25  Hz  and  redundant  traces  editted. 
A  reduction  velocity  of  4.6  km/sec  has  been  used.  The 
topographic  profile  has  a  vertical  exaggeration  of  two 
times.  The  travel  times  for  the  first  arrivals  from  the 
Dinoseis  shots  are  also  plotted. 
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for  the  first  arrivals  quite  well.  Future  interpretation 
efforts  will  start  with  geologic  and  geophysical  cross 
sections  developed  from  borehole  information.  A  computer 
program  based  on  the  psuedospectral  approximation  for  the 
acoustic  and  elastic  wave  equations  (Kosloff  and  Baysal, 
1982  and  Kosloff  et  al.,  1984)  will  be  used  to  model  the 
seismic  wave  field  produced  by  this  experiment.  The 
resulting  cross  section  should  agree  with  the  boreholes, 
gravity  data  and  the  seismic  section,  and  hence  should  be 
an  accurate  representation  of  the  caldera  structure 
across  its  southeast  margin.  The  model  should  predict  the 
seismic  response  of  nuclear  tests  located  near  the 
profile. 
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Appendix 

Derivation  of  the  Fourier  gravity  inversion  solution  with 
density  a  linear  function  of  depth: 

Starting  with  the  expression  for  the  gravitational 
potential  due  to  attraction  from  a  layer  of  material: 


U(r0)  =  G 


c(z)  /  j  r0  -  r |  dV 


where: 


G 

V 

o(z) 

r0 

r 

I 


gravitational  constant 
volume 

linear  density  function 

o(z)  =  c0  +  cz 

position  of  the  observation 

position  within  V 

projection  of  r  onto  the  surface 
z  =  0. 


and  expanding  the  integral  expression, 

„  h(r) 


U(r0)  =  G 


o(z)  /  j rn  -  r|  dz  dS 


D  J 


Taking  the  Fourier  transform  of  the  gravitational 
potential,  the  expression  becomes: 


F[U(r0)]  = 


F[U(r0) ]  =  G 


U(rn)  exp(ik-r)  dS, 


,h(r) 


exp(i|k|  (r0  -  r) )  exp(iKT)  dz  dS( 
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« 

F[U(r0)]  =  2tt  G  exp  ( i  |  k  |  r)  exp(-|k|z0) 

JD 

•  {  Oq  ( exp  ( |  )$  |  h  ( r )  )  -  1)  /  |k|  +  c((|k|h(r)  -  1) 
'  exp(  |  )s|h(r)  )  *  1)  /  ik|2)  dS 


After  Taylor  series  expansion  of  the  exponential  and 
rear  'anging  terms: 


oo 


F[U(r0)  ]  =  2ir  G  exp(-|k|z0)  {aQ  2  \k\n~2  /  n!  F[hn(r)] 

n=l 

oo 

+  c  2  | k| n”3  /  n!  (n  -  1)  F[hn(r)]} 
n=2 

To  include  both  an  upper  and  lower  surface,  this 
expression  becomes: 

00 

F[U(r0)  ]  =  2  IT  G  exp  ( -  |  k  |  z  0 )  {aQ  2  \k\n~2  /  n!  F[hn(r) 

n=l 

oo 

-  gn(r)  ]  +  c  2  1  k |  n“3  /  n I  (n  -  l)  F[hn(j;) 

n=2 

-  gn(r) ] } 


To  obtain  the  expression  for  the  gravity  effect, 
differentiate  both  sides  with  respect  to  z  and  evaluate 


at  Zq  : 


00 


F[Sg(r)]  =  2TT  G  exp(-|k|z0)  { aQ  2  |k|n_1  /  n!  F[hn(r) 

n=l 


-  gn(£)]  +  c  2  | k | n“2  /  n!  (n  -  1)  F[hn(r) 

n=2 

-  gn(£) ] ) 

To  obtain  the  inverse  expression,  solve  for  the  first 
term  of  the  first  infinite  series  as  follows: 


F[  Sg(£)  ]  =2  ir  G  exp(-|Klz0)  o0  F[hn(r)  -  gx,(r)] 


.n, 
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'1 


GO 

+  2n  G  exp(-|k|zQ)  { oQ  Z  |)s|n-1  /  n!  F[hn(r) 

n=2 

00 

-  gn(£)3  +  c  Z  |  k| n-2  /  n!  (n  -  1)  F[hn(r) 

n=2 

-  gn(r) ] ) 

F[hn(r)  -  gn (£)  ]  =  F[Sg(r)]  exp(|fc|z0)  /  (2tr  G  o0) 

co  oo 

-  {  S  |  k| n_1  /  n!  +  (c  /  aQ)  Z  |  k| n-2  /  n! 

n-2  n=2 


•  (n  -  1) }  F[hn (r)  -  gn(r) ] 
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